home *** CD-ROM | disk | FTP | other *** search
/ Skunkware 5 / Skunkware 5.iso / lib / calc / ellip.cal < prev    next >
Text File  |  1995-07-17  |  5KB  |  173 lines
  1. /*
  2.  * Copyright (c) 1993 David I. Bell
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * Attempt to factor numbers using elliptic functions.
  7.  *     y^2 = x^3 + a*x + b   (mod N).
  8.  *
  9.  * Many points (x,y) (mod N) are found that solve the above equation,
  10.  * starting from a trivial solution and 'multiplying' that point together
  11.  * to generate high powers of the point, looking for such a point whose
  12.  * order contains a common factor with N.  The order of the group of points
  13.  * varies almost randomly within a certain interval for each choice of a
  14.  * and b, and thus each choice provides an independent opportunity to
  15.  * factor N.  To generate a trivial solution, a is chosen and then b is
  16.  * selected so that (1,1) is a solution.  The multiplication is done using
  17.  * the basic fact that the equation is a cubic, and so if a line hits the
  18.  * curve in two rational points, then the third intersection point must
  19.  * also be rational.  Thus by drawing lines between known rational points
  20.  * the number of rational solutions can be made very large.  When modular
  21.  * arithmetic is used, solving for the third point requires the taking of a
  22.  * modular inverse (instead of division), and if this fails, then the GCD
  23.  * of the failing value and N provides a factor of N.  This description is
  24.  * only an approximation, read "A Course in Number Theory and Cryptography"
  25.  * by Neal Koblitz for a good explanation.
  26.  *
  27.  * factor(iN, ia, B, force)
  28.  *    iN is the number to be factored.
  29.  *    ia is the initial value of a in the equation, and each successive
  30.  *    value of a is an independent attempt at factoring (default 1).
  31.  *    B is the limit of the primes that make up the high power that the
  32.  *    point is raised to for each factoring attempt (default 100).
  33.  *    force is a flag to attempt to factor numbers even if they are
  34.  *    thought to already be prime (default FALSE).
  35.  *
  36.  * Making B larger makes the power the point being raised to contain more
  37.  * prime factors, thus increasing the chance that the order of the point
  38.  * will be made up of those factors.  The higher B is then, the greater
  39.  * the chance that any individual attempt will find a factor.  However,
  40.  * a higher B also slows down the number of independent functions being
  41.  * examined.  The order of the point for any particular function might
  42.  * contain a large prime and so won't succeed even for a really large B,
  43.  * whereas the next function might have an order which is quickly found.
  44.  * So you want to trade off the depth of a particular search with the
  45.  * number of searches made.  For example, for factoring 30 digits, I make
  46.  * B be about 1000 (probably still too small).
  47.  *
  48.  * If you have lots of machines available, then you can run parallel
  49.  * factoring attempts for the same number by giving different starting
  50.  * values of ia for each machine (e.g. 1000, 2000, 3000).
  51.  *
  52.  * The output as the function is running is (occasionally) the value of a
  53.  * when a new function is started, the prime that is being included in the
  54.  * high power being calculated, and the current point which is the result
  55.  * of the powers so far.
  56.  *
  57.  * If a factor is found, it is returned and is also saved in the global
  58.  * variable f.  The number being factored is also saved in the global
  59.  * variable N.
  60.  */
  61.  
  62. obj point {x, y};
  63. global    N;        /* number to factor */
  64. global    a;        /* first coefficient */
  65. global    b;        /* second coefficient */
  66. global    f;        /* found factor */
  67.  
  68.  
  69. define factor(iN, ia, B, force)
  70. {
  71.     local    C, x, p;
  72.  
  73.     if (!force && ptest(iN, 50))
  74.         return 1;
  75.     if (isnull(B))
  76.         B = 100;
  77.     if (isnull(ia))
  78.         ia = 1;
  79.     obj point x;
  80.     a = ia;
  81.     b = -ia;
  82.     N = iN;
  83.     C = isqrt(N);
  84.     C = 2 * C + 2 * isqrt(C) + 1;
  85.     f = 0;
  86.     while (f == 0) {
  87.         print "A =", a;
  88.         x.x = 1;
  89.         x.y = 1;
  90.         print 2, x;
  91.         x = x ^ (2 ^ (highbit(C) + 1));
  92.         for (p = 3; ((p < B) && (f == 0)); p += 2) {
  93.             if (!ptest(p, 1))
  94.                 continue;
  95.             print p, x;
  96.             x = x ^ (p ^ ((highbit(C) // highbit(p)) + 1));
  97.         }
  98.         a++;
  99.         b--;
  100.     }
  101.     return f;
  102. }
  103.  
  104.  
  105. define point_print(p)
  106. {
  107.     print "(" : p.x : "," : p.y : ")" :;
  108. }
  109.  
  110.  
  111. define point_mul(p1, p2)
  112. {
  113.     local    r, m;
  114.  
  115.     if (p2 == 1)
  116.         return p1;
  117.     if (p1 == p2)
  118.         return point_square(&p1);
  119.     obj point r;
  120.     m = (minv(p2.x - p1.x, N) * (p2.y - p1.y)) % N;
  121.     if (m == 0) {
  122.         if (f == 0)
  123.             f = gcd(p2.x - p1.x, N);
  124.         r.x = 1;
  125.         r.y = 1;
  126.         return r;        
  127.     }
  128.     r.x = (m^2 - p1.x - p2.x) % N;
  129.     r.y = ((m * (p1.x - r.x)) - p1.y) % N;
  130.     return r;
  131. }
  132.  
  133.  
  134. define point_square(p)
  135. {
  136.     local    r, m;
  137.  
  138.     obj point r;
  139.     m = ((3 * p.x^2 + a) * minv(p.y << 1, N)) % N;
  140.     if (m == 0) {
  141.         if (f == 0)
  142.             f = gcd(p.y << 1, N);
  143.         r.x = 1;
  144.         r.y = 1;
  145.         return r;
  146.     }
  147.     r.x = (m^2 - p.x - p.x) % N;
  148.     r.y = ((m * (p.x - r.x)) - p.y) % N;
  149.     return r;
  150. }
  151.  
  152.  
  153. define point_pow(p, pow)
  154. {
  155.     local bit, r, t;
  156.  
  157.     r = 1;
  158.     if (isodd(pow))
  159.         r = p;
  160.     t = p;
  161.     for (bit = 2; ((bit <= pow) && (f == 0)); bit <<= 1) {
  162.         t = point_square(&t);
  163.         if (bit & pow)
  164.             r = point_mul(&t, &r);
  165.     }
  166.     return r;
  167. }
  168.  
  169. global lib_debug;
  170. if (lib_debug >= 0) {
  171.     print "factor(N, I, B, force) defined";
  172. }
  173.